%     VIRTUAL LAMINATE - программа для получения эффективных свойств слоистых пластин
%     Copyright (C) 2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function Write_Input_Thermal(Job_Name, Work_Directory, NODES, ELEMENTS, RVE_Coordinates, Node_Merge_Tolerance, Ply_Elements, Materials, Stacking, Element_Order)

global ErrCode;

node_coord_format = '%.12e';
orientation_format = '%.12e';
mat_prop_format = '%.12e';
equation_format = '%.12e';
bc_format = '%.12e';

Part_Name = 'LAMINATE_RVE';
Model_Name = Job_Name;

% Проверка нумерации узлов
for i=1:size(NODES,1)
    if (NODES(i,1) ~= i)
        fprintf (1, '\nОшибка функции Write_Input_Thermal: некорректная нумерация узлов!\n');
        ErrCode = -1;
        return;
    end
end

% Проверка нумерации элементов
for i=1:size(ELEMENTS,1)
    if (ELEMENTS(i,1) ~= i)
        fprintf (1, '\nОшибка функции Write_Input_Thermal: некорректная нумерация элементов!\n');
        ErrCode = -1;
        return;
    end
end

File_name = Job_Name;
File_path = [Work_Directory '\' File_name '.inp'];

Fout = fopen(File_path, 'wt');

if (Fout < 0)
    fprintf (1, '\nОшибка функции Write_Input_Thermal: не удалось открыть для записи файл [%s] !\n', File_path);
    ErrCode = -1;
    return;
end

fprintf (Fout, '*Heading\n');

fprintf (Fout, '** Job name: %s Model name: %s\n', Job_Name, Model_Name);
fprintf (Fout, '** Generated by: Virtual_Laminate\n');

fprintf (Fout, '*Preprint, echo=NO, model=NO, history=NO, contact=NO\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** PARTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Part, name=%s\n', Part_Name);

fprintf (Fout, '*Node\n');

for i=1:size(NODES, 1)
    fprintf (Fout, ['%i, ' node_coord_format ', ' node_coord_format ', ' node_coord_format '\n'], NODES(i,1), NODES(i,2), NODES(i,3), NODES(i,4));
end

switch Element_Order
    
    case 'LINEAR'
        
        fprintf (Fout, '*Element, type=DC3D8\n');

        for i=1:size(ELEMENTS,1)
            fprintf (Fout, '%i, %i, %i, %i, %i, %i, %i, %i, %i\n', ELEMENTS(i,1), ELEMENTS(i,5), ELEMENTS(i,6), ELEMENTS(i,7), ...
                                                                   ELEMENTS(i,8), ELEMENTS(i,9), ELEMENTS(i,10), ELEMENTS(i,11), ...
                                                                   ELEMENTS(i,12));
        end
        
    case 'QUADRATIC'
        
        fprintf (Fout, '*Element, type=DC3D20\n');

        for i=1:size(ELEMENTS,1)
            fprintf (Fout, '%i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i, %i,\n', ELEMENTS(i,1), ELEMENTS(i,5), ELEMENTS(i,6), ELEMENTS(i,7), ...
                                                                                               ELEMENTS(i,8), ELEMENTS(i,9), ELEMENTS(i,10), ELEMENTS(i,11), ...
                                                                                               ELEMENTS(i,12), ELEMENTS(i,13), ELEMENTS(i,14), ELEMENTS(i,15), ...
                                                                                               ELEMENTS(i,16), ELEMENTS(i,17), ELEMENTS(i,18), ELEMENTS(i,19));
            fprintf (Fout, '   %i, %i, %i, %i, %i\n', ELEMENTS(i,20), ELEMENTS(i,21), ELEMENTS(i,22), ELEMENTS(i,23), ELEMENTS(i,24));
        end
    
end

fprintf (Fout, '*Nset, nset=ALL_NODES, generate\n');
fprintf (Fout, ' %i, %i, 1\n', 1, NODES(size(NODES,1),1));

fprintf (Fout, '*Elset, elset=ALL_ELEMENTS, generate\n');
fprintf (Fout, ' %i, %i, 1\n', 1, ELEMENTS(size(ELEMENTS,1),1));

N_digits = numel(num2str(Stacking(size(Stacking,1),1)));

format_line = ['*Elset, elset=Ply_%.' sprintf('%i', N_digits) 'i\n'];

N_previous_items = 0;

for i=1:size(Stacking,1)
        
    fprintf (Fout, format_line, Stacking(i,1));
    
    Write_Set(Fout, ELEMENTS(N_previous_items + 1 : N_previous_items + Ply_Elements(i), 1));
    
    N_previous_items = N_previous_items + Ply_Elements(i);
    
end

format_line = ['*Orientation, name=Ori_Coord_%.' sprintf('%i', N_digits) 'i\n'];

format_line_2 = [orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format ', ' orientation_format '\n'];

Origin_High = 0;

for i=1:size(Stacking,1)
    
    fprintf (Fout, format_line, Stacking(i,1));
    
    fprintf (Fout, format_line_2,  cosd(Stacking(i,2)), sind(Stacking(i,2)), Origin_High + Materials{Stacking(i,3)}.t/2, ...
                                  -sind(Stacking(i,2)), cosd(Stacking(i,2)), Origin_High + Materials{Stacking(i,3)}.t/2, ...
                                                                       0, 0, Origin_High + Materials{Stacking(i,3)}.t/2);
    fprintf (Fout, '3, 0.\n');
    
    Origin_High = Origin_High + Materials{Stacking(i,3)}.t;
    
end

format_line = ['** Section: Section_Ply_%.' sprintf('%i', N_digits) 'i\n'];

format_line_2 = ['*Solid Section, elset=Ply_%.' sprintf('%i', N_digits) 'i, orientation=Ori_Coord_%.' sprintf('%i', N_digits) 'i, material=Material_%i\n'];

for i=1:size(Stacking,1)
    
    fprintf (Fout, format_line, Stacking(i,1));
    
    fprintf (Fout, format_line_2, Stacking(i,1), Stacking(i,1), Stacking(i,3));
    
    fprintf (Fout, ',\n');
    
end

fprintf (Fout, '*End Part\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** ASSEMBLY\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Assembly, name=Assembly\n');
fprintf (Fout, '**\n');
fprintf (Fout, '*Instance, name=%s-1, part=%s\n', Part_Name, Part_Name);
fprintf (Fout, '*End Instance\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node\n');
fprintf (Fout, '      1,           0.,           0.,           0.\n');
fprintf (Fout, '      2,           0.,           0.,           0.\n');
fprintf (Fout, '      3,           0.,           0.,           0.\n');

fprintf (Fout, '*Nset, nset=RP-ThermGrad_x\n');
fprintf (Fout, ' 1,\n');

fprintf (Fout, '*Nset, nset=RP-ThermGrad_y\n');
fprintf (Fout, ' 2,\n');

fprintf (Fout, '*Nset, nset=RP-ThermGrad_z\n');
fprintf (Fout, ' 3,\n');

RVE_AG_NegX_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,1)+Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);
RVE_AG_PosX_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(1,2)+Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);

RVE_AG_NegY_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,1)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);
RVE_AG_PosY_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,2)-Node_Merge_Tolerance, RVE_Coordinates(2,2)+Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);

RVE_AG_NegZ_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)-Node_Merge_Tolerance, RVE_Coordinates(3,1)-Node_Merge_Tolerance, RVE_Coordinates(3,1)+Node_Merge_Tolerance);
RVE_AG_PosZ_RVE_AGSORT = Generate_Nset(NODES, RVE_Coordinates(1,1)-Node_Merge_Tolerance, RVE_Coordinates(1,2)-Node_Merge_Tolerance, RVE_Coordinates(2,1)-Node_Merge_Tolerance, RVE_Coordinates(2,2)-Node_Merge_Tolerance, RVE_Coordinates(3,2)-Node_Merge_Tolerance, RVE_Coordinates(3,2)+Node_Merge_Tolerance);

if (isempty(RVE_AG_NegX_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_NegX_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_NegY_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_NegY_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_NegZ_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_NegZ_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_PosX_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_PosX_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_PosY_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_PosY_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (isempty(RVE_AG_PosZ_RVE_AGSORT))
    fprintf (1, '\nОшибка: не удалось найти узлы для сэта RVE_AG_PosZ_RVE_AGSORT !\n');
    fprintf (1, 'Измените параметры построения сетки.\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (numel(RVE_AG_NegX_RVE_AGSORT) ~= numel(RVE_AG_PosX_RVE_AGSORT))
    fprintf (1, '\nОшибка: размер сэта RVE_AG_NegX_RVE_AGSORT не совпадает с размером сэта RVE_AG_PosX_RVE_AGSORT !\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (numel(RVE_AG_NegY_RVE_AGSORT) ~= numel(RVE_AG_PosY_RVE_AGSORT))
    fprintf (1, '\nОшибка: размер сэта RVE_AG_NegY_RVE_AGSORT не совпадает с размером сэта RVE_AG_PosY_RVE_AGSORT !\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

if (numel(RVE_AG_NegZ_RVE_AGSORT) ~= numel(RVE_AG_PosZ_RVE_AGSORT))
    fprintf (1, '\nОшибка: размер сэта RVE_AG_NegZ_RVE_AGSORT не совпадает с размером сэта RVE_AG_PosZ_RVE_AGSORT !\n');
    fclose(Fout);
    ErrCode = -1;
    return;
end

RVE_AG_PosX_RVE_AGSORT = Match_Nodes(RVE_AG_PosX_RVE_AGSORT, RVE_AG_NegX_RVE_AGSORT, NODES, 1, Node_Merge_Tolerance);

if (ErrCode < 0)
    fclose(Fout);
    return;
end

RVE_AG_PosY_RVE_AGSORT = Match_Nodes(RVE_AG_PosY_RVE_AGSORT, RVE_AG_NegY_RVE_AGSORT, NODES, 2, Node_Merge_Tolerance);

if (ErrCode < 0)
    fclose(Fout);
    return;
end

RVE_AG_PosZ_RVE_AGSORT = Match_Nodes(RVE_AG_PosZ_RVE_AGSORT, RVE_AG_NegZ_RVE_AGSORT, NODES, 3, Node_Merge_Tolerance);

if (ErrCode < 0)
    fclose(Fout);
    return;
end

fprintf (Fout, '*NSET, NSET=RVE_AG_NegX_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_NegX_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_PosX_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_PosX_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_NegY_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_NegY_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_PosY_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_PosY_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_NegZ_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_NegZ_RVE_AGSORT);

fprintf (Fout, '*NSET, NSET=RVE_AG_PosZ_RVE_AGSORT, UNSORTED, instance=%s-1\n', Part_Name);
Write_Set(Fout, RVE_AG_PosZ_RVE_AGSORT);

fprintf (Fout, '**\n');
fprintf (Fout, '** EQUATIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Equation\n');

fprintf (Fout, '3\n');
fprintf (Fout, ['RVE_AG_PosX_RVE_AGSORT, 11, -1.000000e+00, RVE_AG_NegX_RVE_AGSORT, 11, 1.000000e+00, RP-ThermGrad_x, 11, ' equation_format '\n'], RVE_Coordinates(1,2) - RVE_Coordinates(1,1));
fprintf (Fout, '3\n');
fprintf (Fout, ['RVE_AG_PosY_RVE_AGSORT, 11, -1.000000e+00, RVE_AG_NegY_RVE_AGSORT, 11, 1.000000e+00, RP-ThermGrad_y, 11, ' equation_format '\n'], RVE_Coordinates(2,2) - RVE_Coordinates(2,1));
fprintf (Fout, '3\n');
fprintf (Fout, ['RVE_AG_PosZ_RVE_AGSORT, 11, -1.000000e+00, RVE_AG_NegZ_RVE_AGSORT, 11, 1.000000e+00, RP-ThermGrad_z, 11, ' equation_format '\n'], RVE_Coordinates(3,2) - RVE_Coordinates(3,1));

fprintf (Fout, '*End Assembly\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** MATERIALS\n');
fprintf (Fout, '**\n');

Mat_Numbers = unique(Stacking(:,3));

for i=1:numel(Mat_Numbers)
    
    fprintf (Fout, '*Material, name=Material_%i\n', Mat_Numbers(i));
    fprintf (Fout, '*Conductivity, type=ORTHO\n');
    fprintf (Fout, [mat_prop_format ', ' mat_prop_format ', ' mat_prop_format '\n'], Materials{Mat_Numbers(i)}.lambda11, Materials{Mat_Numbers(i)}.lambda22, Materials{Mat_Numbers(i)}.lambda33);
    fprintf (Fout, '*Density\n');
    fprintf (Fout, [mat_prop_format ',\n'], Materials{Mat_Numbers(i)}.rho);
    fprintf (Fout, '*Specific Heat\n');
    fprintf (Fout, [mat_prop_format ',\n'], Materials{Mat_Numbers(i)}.Cp);

end

fprintf (Fout, '** ----------------------------------------------------------------\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** STEP: HomogK_X-Initial\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Step, name=HomogK_X-Initial, nlgeom=NO\n');
fprintf (Fout, 'RP_grad\n');
fprintf (Fout, '*Heat Transfer, steady state, deltmx=0.\n');
fprintf (Fout, '1e-08, 1e-08, 1e-13, 1e-08, \n');

fprintf (Fout, '**\n');
fprintf (Fout, '** BOUNDARY CONDITIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary\n');
fprintf (Fout, ['RP-ThermGrad_x, 11, 11, ' bc_format '\n'], 1/(RVE_Coordinates(1,2) - RVE_Coordinates(1,1)));

fprintf (Fout, '** Name: HomogK_X-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary\n');
fprintf (Fout, 'RP-ThermGrad_y, 11, 11\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary\n');
fprintf (Fout, 'RP-ThermGrad_z, 11, 11\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** OUTPUT REQUESTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Restart, write, frequency=0\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: ThermalFieldOutput\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field\n');
fprintf (Fout, '*Node Output\n');
fprintf (Fout, 'NT, RFL\n');

fprintf (Fout, '*Element Output, directions=YES\n');
fprintf (Fout, 'FLUXS, HFL, IVOL, TEMP\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: F-Output-1\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field, variable=PRESELECT\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_x_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, history\n');
fprintf (Fout, '*Node Output, nset=RP-ThermGrad_x\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_y_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_y\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_z_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_z\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '*EL PRINT, ELSET=%s-1.ALL_ELEMENTS, POSITION=INTEGRATION POINTS, FREQUENCY=99999\n', Part_Name);
fprintf (Fout, 'HFL, IVOL\n');

fprintf (Fout, '*End Step\n');

fprintf (Fout, '** ----------------------------------------------------------------\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** STEP: HomogK_Y-Initial\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Step, name=HomogK_Y-Initial, nlgeom=NO\n');
fprintf (Fout, 'RP_grad\n');
fprintf (Fout, '*Heat Transfer, steady state, deltmx=0.\n');
fprintf (Fout, '1e-08, 1e-08, 1e-13, 1e-08, \n');

fprintf (Fout, '**\n');
fprintf (Fout, '** BOUNDARY CONDITIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_X-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_x, 11, 11\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, ['RP-ThermGrad_y, 11, 11, ' bc_format '\n'], 1/(RVE_Coordinates(2,2) - RVE_Coordinates(2,1)));

fprintf (Fout, '** Name: HomogK_Y-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_z, 11, 11\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** OUTPUT REQUESTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Restart, write, frequency=0\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: ThermalFieldOutput\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field\n');
fprintf (Fout, '*Node Output\n');
fprintf (Fout, 'NT, RFL\n');

fprintf (Fout, '*Element Output, directions=YES\n');
fprintf (Fout, 'FLUXS, HFL, IVOL, TEMP\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: F-Output-1\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field, variable=PRESELECT\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_x_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, history\n');
fprintf (Fout, '*Node Output, nset=RP-ThermGrad_x\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_y_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_y\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_z_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_z\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '*EL PRINT, ELSET=%s-1.ALL_ELEMENTS, POSITION=INTEGRATION POINTS, FREQUENCY=99999\n', Part_Name);
fprintf (Fout, 'HFL, IVOL\n');

fprintf (Fout, '*End Step\n');

fprintf (Fout, '** ----------------------------------------------------------------\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** STEP: HomogK_Z-Initial\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Step, name=HomogK_Z-Initial, nlgeom=NO\n');
fprintf (Fout, 'RP_grad\n');
fprintf (Fout, '*Heat Transfer, steady state, deltmx=0.\n');
fprintf (Fout, '1e-08, 1e-08, 1e-13, 1e-08, \n');

fprintf (Fout, '**\n');
fprintf (Fout, '** BOUNDARY CONDITIONS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Y-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');

fprintf (Fout, '** Name: HomogK_Z-Initial-dTX Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_x, 11, 11\n');

fprintf (Fout, '** Name: HomogK_Z-Initial-dTY Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, 'RP-ThermGrad_y, 11, 11\n');

fprintf (Fout, '** Name: HomogK_Z-Initial-dTZ Type: Temperature\n');
fprintf (Fout, '*Boundary, op=NEW\n');
fprintf (Fout, ['RP-ThermGrad_z, 11, 11, ' bc_format '\n'], 1/(RVE_Coordinates(3,2) - RVE_Coordinates(3,1)));

fprintf (Fout, '**\n');
fprintf (Fout, '** OUTPUT REQUESTS\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Restart, write, frequency=0\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: ThermalFieldOutput\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field\n');
fprintf (Fout, '*Node Output\n');
fprintf (Fout, 'NT, RFL\n');

fprintf (Fout, '*Element Output, directions=YES\n');
fprintf (Fout, 'FLUXS, HFL, IVOL, TEMP\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** FIELD OUTPUT: F-Output-1\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, field, variable=PRESELECT\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_x_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Output, history\n');
fprintf (Fout, '*Node Output, nset=RP-ThermGrad_x\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_y_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_y\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '**\n');
fprintf (Fout, '** HISTORY OUTPUT: RP-ThermGrad_z_History\n');
fprintf (Fout, '**\n');

fprintf (Fout, '*Node Output, nset=RP-ThermGrad_z\n');
fprintf (Fout, 'CFL, NT, RFL\n');

fprintf (Fout, '*EL PRINT, ELSET=%s-1.ALL_ELEMENTS, POSITION=INTEGRATION POINTS, FREQUENCY=99999\n', Part_Name);
fprintf (Fout, 'HFL, IVOL\n');

fprintf (Fout, '*End Step\n');

fclose(Fout);
